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In  this  project,  the  uniquely-developed  TOhoku  Mixed  Basis  Orbital  ab  initio  calculation  package 
(TOMBO)  was  used  in  simulations  of  electron  and  thermal  transport  in  a  small  molecule  organic-based 
circuit  and  in  a  GeSe  supercell,  respectively.  For  the  small  molecule,  application  of  a  bias  lead  to  a 
non-equilibrium  state  for  which  Density  Functional  Theory  (DFT)  is  not  directly  applicable.  The 
simulation  proved  useful  by  modifying  TOMBO  by  separating  the  whole  system  (circuit)  into  three  parts: 
right  and  left  semi-infinitive  leads,  and  finite  organic  molecule.  Since  TOMBO  uses  two  different  types  of 
basis  sets  (plane-waves  and  atomic  orbitals)  maximally  localized  Wannier  function  (Wannier90)  was 
introduced  to  realize  the  separation.  Since  the  two  semi-infinitive  leads  are  very  close  to  their  equilibrium 
state,  non-equilibrium  Green  function  (NEGF)  was  used  to  represent  them  by  self-energy  terms. 
Introduction  of  bias  affected  periodic  boundary  conditions,  but  were  satisfied  by  introducing  a  potential 
drop  in  the  vacuum.  Simulation  of  thermal  conductivity  versus  temperature  was  considered  for  a  6x6x6 
k-mesh  involving  a  crystalline,  orthorhombic  GeSe  supercell.  Its  thermal  conductivity  tensor  is  diagonal, 
with  three  different  values  along  the  3  principal  axes.  This  work  examined  the  third  trace  of  that  tensor 
and  the  contribution  of  the  different  branches.  Acoustic  branches  usually  contribute  most  to  the  thermal 
conductivity  as  their  group  velocity  and  lifetimes  are  largest.  In  this  case,  each  of  the  acoustic  modes 
contributes  to  about  20%  of  the  total  value.  The  first  3  branches  contribute  to  almost  half  the  total  thermal 
conductivity. 
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1.  First-principles  calculation  of  Electron  transport  through  organic 
molecule 


Y.Y.  Liang,  Y.  Kawazoe 

New  Industry  Creation  Hatchery  Center,  Tohoku  University,  Aobaku,  Sendai,  980-8579  Japan 


Motivation: 


The  motivation  of  our  study  is  to  realize  the  organic  molecular  conductivity 
computation  based  on  the  first-principles  calculation.  Organic  molecules  are 
considered  as  one  of  the  prospective  candidates  of  the  next-generation  circuit.  As 
shown  in  the  figure,  the  molecule  connects  to  the  two  leads.  If  the  bias  is  applied, 
the  electron  current  flows  through  the  molecule.  The  I-V  curve  is  calculated  self- 
consistently.  The  computation  is  realized  by  self-developed  package  TOMBO. 


I 


Fig.  1.  Benzene  is  connected  with  two  Au  leads 


Different  from  the  traditional  density  functional  theory  [DFT)  method,  the 
present  system  is  divided  into  two  infinitive  leads  connecting  with  organic 
molecules  [shown  in  Fig.  1.).  The  whole  is  in  non-equilibrium  after  applying  the 
bias.  Therefore,  we  cannot  apply  the  DFT  calculations  directly.  To  achieve  the 
computation  by  using  TOMBO,  some  modifications  are  required. 

Methods: 


Original  Algorithm  in  TOMBO  New  algorithm  in  TOMBO 


Fig.  2.  The  algorithms  used  in  traditional  DFT  calculations  and  transport 

calculations. 

In  traditional  DFT  calculations,  the  equations  are  solved  self-consistently  [shown 
in  Fig.  2.).  Firstly,  based  on  DFT,  the  initial  density  distribution  results  into  the 
Hamiltonian  of  the  system.  Then,  from  statistical  mechanics,  diagonalization  of 
the  Hamiltonian  achieves  new  distribution.  Generally  speaking,  this  new 


distribution  is  different  from  the  previous  one.  Applying  the  self-consistent 
process,  the  system  is  calculated  until  the  criteria  are  satisfied. 


However,  in  present  case,  since  the  introduction  of  two  semi-infinitive  leads,  the 
whole  system,  which  has  no  periodicity,  is  infinitive.  The  traditional  DFT  fails  to 
deal  with  the  infinitive  system.  Therefore,  non-equilibrium  Green  function 
(NEGF)  is  introduced  and  two  semi-infinitive  leads  are  represented  by  self¬ 
energy  terms. 

Compared  with  the  traditional  DFT  process,  the  first  step  is  same:  the 
Hamiltonian  comes  from  the  density  distribution.  But  the  second  step  is  replaced, 
and  the  new  density  distribution  stems  from  the  NEGF  rather  than  the  statistical 
mechanics. 

NEGF: 

The  Hamiltonian  (H]  of  the  system  reads, 


H  =  Y  (kicliCki  +  Y  Etd\dl  +  Y  Vki,lckidl  +  c.c 


The  first  term  is  the  energy  levels  of  the  leads,  the  second  term  is  the  energy 
levels  of  the  molecules  connecting  with  the  leads  and  the  last  two  term  stand  for 
the  coupling  between  the  molecules  and  the  leads. 

The  self-energy  term  is, 


where  giR  is  the  surface  Green  function.  After  achievement  of  self-energy,  we  can 
get  the  new  density  matrix, 


The  density  matrix  and  the  Hamiltonian  are  solved  self-consistently,  and  the 
current  is  calculated  if  the  criteria  are  satisfied. 


I  =  ~h  J  dET1(rlGRr.2oA)(ME)  -  me)) 


where  ft  is  the  Fermi  distribution  of  the  ith  lead  and 


which  represents  the  broaden  of  the  energy  levels. 


Wannier  Function: 


The  basis  sets  used  in  TOMBO  is  mixed  basis  sets. 


tGr 


In  the  mixed  basis  sets,  there  are  two  parts.  One  is  plane  waves  basis  sets,  the 
other  is  atomic  basis  sets. 

The  plane  waves  are  non-localized.  To  divide  the  whole  system  into  three  parts, 
we  have  to  introduce  Wannier  function.  In  our  program,  we  will  use  maximally 
localized  Wannier  function  (MLWF). 

If  Wannier  function  is  achieved,  we  can  express  the  Hamiltonian  operator  in  real 
space  by  using  the  definition: 

'Rn)-^ffiu^y"Sdk 

where  Uk  is  the  unitary  matrix.  The  strategy  consists  of  choosing  the  U  that 
minimize  the  sum  of  the  following  spreads  of  the  Wannier  functions  about  their 
centers: 


Q  =  ^  (On  |  r2 1  On^  —  (On  \  r  |  0w)2J  =  ^  —  rn  1 

71  71  L  J 

The  matrix  elements  of  the  position  operator  between  Wannier  functions  may  be 
expressed  in  reciprocal  space  as: 

(Rn\r\0m)  =  i  f  e*'R  (unk  |  V*  |  umk  )dk 


(2tt) 


(Rti  |  r2  1 0 m'j  = 


(2rrj 


fe*'X 


If  the  Wannier  function  is  achieved,  we  can  get  the  self-energy  and  calculate  the 
current  through  the  molecule. 


Outlook: 


Since  the  periodic  boundary  condition  should  be  valid  in  the  calculation,  the 
application  of  the  electric  field  breaks  the  translation  invariance.  To  recover  the 
periodic  boundary  condition,  the  saw-tooth  like  potential  should  be  applied  in 
our  calculation.  So  far,  the  application  of  the  electric  field  is  not  available  in 
TOMBO,  we  have  to  modify  the  source  code  and  make  it  possible. 

Conclusion: 


1]  The  whole  system  should  be  separated  into  three  parts,  left  semi¬ 
infinitive  lead,  finite  organic  molecule,  and  the  right  semi-infinitive  lead. 
Since  two  different  types  of  basis  sets,  namely  plane-waves  and  atomic 
orbitals,  are  used  in  TOMBO,  we  introduce  the  maximally  localized 
Wannier  function  (Wannier90),  to  realize  the  separation. 

2)  The  infinitive  system  cannot  be  calculated  directly.  Given  the  fact  that  two 
leads  are  very  close  to  their  equilibrium  state,  non-equilibrium  Green 
function  (NEGF)  is  used,  and  two  semi-infinitive  leads  are  represented  by 
self-energy  terms.  These  self-energy  terms  have  the  same  dimensions  as 


the  organic  molecule.  Therefore,  the  infinitive  system  is  equivalent  to  the 
finite  organic  molecule  plus  two  finite  self-energy  terms. 

3)  In  the  calculation,  the  periodical  boundary  condition  should  be  satisfied. 
However,  if  the  bias  is  applied,  the  potential  drop  is  found  in  the  two  leads 
to  counter-balance  the  applied  bias.  Unfortunately,  this  is  not  the  real  case. 
To  solve  this  problem,  we  introduce  the  potential  drop  in  the  vacuum. 
Therefore,  electrostatic  potential  of  two  leads  is  homogeneous,  and  the 
leads  can  be  replaced  by  two  self-energy  terms. 


1.  First-principles  calculation  of  thermal  transport  in  GeSe 


Keivan  Esfarjani1,  Yunye  Liang2  and  Yoshiyuki  Kawazoe2 

1  Department  of  Mechanical  Engineering,  Rutgers  University, 
Piscataway,  NJ,  USA 

2  Institute  for  Materials  Research,  Tohoku  University,  Sendai,  Japan 

GeSe  as  a  IV-VI  chalcogenide  material,  has  good  potential  as  a 
thermoelectric  material.  Our  goal  in  this  project  is  to  apply  our 
newly-developed  methodology[l]  to  calculate  its  phonon  dispersion 
and  thermal  conductivity. 

In  the  previous  reports,  from  the  first-principles  DFT  calculations  in  a 
large  supercell  containing  480  atoms,  we  obtained  the  phonon 
dispersion  and  showed  that  a  large  number  of  neighbor  interactions 
needs  to  be  included  in  order  to  produce  the  correct  phonon 
dispersion. 

Due  to  the  low  symmetry  of  the  structure  and  the  nature  of  the 
chalcogenide  element,  the  phonon  calculations  did  not  converge  until 
the  range  of  interactions  were  increased  up  to  8  A0 .  This  includes  up 
to  41 st  nearest  neighbors  in  the  harmonic  force  constants.  Cubic  FCs 
were  included  up  to  5th  shell  and  quartic  up  to  the  second  shell. 

Figure  1  shows  the  decay  of  harmonic  spring  constants  versus 
distance  when  up  to  76  neighbor  shells,  corresponding  to  a  pair 
distance  of  9.42  A0  are  included  in  the  fitting.  This  figure  shows  that 
changing  the  number  of  included  force  constants  does  not  change 
their  actually  fitted  value.  Due  to  the  heavy  size  of  the  calculations, 
the  k-point  summation  in  the  calculation  of  the  phonon  lifetimes 
included  a  6x6x6  mesh  in  the  full  Brillouin  zone.  Results  presented  in 
this  report  refer  to  this  mesh.  We  need  however  to  perform 
calculations  with  a  larger  k-mesh  in  order  to  obtain  results  which  are 
converged  with  respect  to  the  number  of  kpoints.  Usually  an 
extrapolation  procedure[l]  needs  to  be  performed  in  order  to  reach 
the  converged  value. 
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Figure  2:  Total  and  partial  thermal  conductivity  versus 
temperature.  The  contribution  of  the  3  acoustic  modes  and 
their  sum  is  being  compared  to  the  total  thermal 
conductivity. 


bic,  the  thermal  conductivity  tensor  is  diagonal,  with  three  different 
values  along  the  3  principal  axes.  Here  we  are  reporting  the  third  of 
the  trace  of  that  tensor.  The  contribution  of  different  branches  can 
also  be  considered.  Usually  the  acoustic  branches  contribute  most  to 
the  thermal  conductivity  as  their  group  velocity  as  well  as  their 
lifetimes  are  largest.  In  this  case,  each  of  the  acoustic  modes 
contributes  to  about  20%  of  the  total  value.  As  can  be  seen  from  this 
figure,  the  first  3  branches  contribute  to  almost  half  the  total  thermal 
conductivity. 


In  the  future,  we  will  examine  the  convergence  of  the  results  with 
respect  to  the  fineness  of  the  k-mesh  chosen  in  the  first  Brillouin 
zone.  These  are  very  heavy  calculations  which  scale  as  the  cube  of  the 
number  of  modes  (=24  in  this  case).  Large  unit  cells  would  require 
much  more  computational  resources  than  the  more  simple  primitive 
cells  which  contain  one  or  two  atoms.  Once  the  convergence  is 
achieved,  phonon  relaxation  times  and  mean  free  paths  will  also  be 
calculated  and  displayed.  Their  contribution  to  the  heat  transport 
will  accordingly  be  discussed. 
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